function spikes=izhikevich_fs(stim_current)

C=20;vr=-55;vt=-40;k=1;
a=0.15;b=8;c=-55;d=200;
vpeak=25;
T=1000;tau=1;


spikes = 0;
n=round(T/tau);
v=vr*ones(1,n); u=0*v;
I=[stim_current*ones(1,n)];

for i = 1:n-1
    v(i+1)=v(i)+tau*(k*(v(i)-vr)*(v(i)-vt)-u(i)+I(i))/C;
    u(i+1)=u(i)+tau*a*(b*(v(i)-vr)-u(i));
    if v(i+1)>=vpeak
        v(i)=vpeak;
        v(i+1)=c;
        u(i+1)=u(i+1)+d;
        spikes = spikes + 1;
    end
end
plot(tau*(1:n),[v' u']);
drawnow
